clip-02-05

clip-02-05.jl

Load Julia packages (libraries) needed for the snippets in chapter 0

using StatisticalRethinking, Optim
gr(size=(600,300));
Plots.GRBackend()

snippet 3.2

p_grid = range(0, step=0.001, stop=1)
prior = ones(length(p_grid))
likelihood = [pdf(Binomial(9, p), 6) for p in p_grid]
posterior = likelihood .* prior
posterior = posterior / sum(posterior)
samples = sample(p_grid, Weights(posterior), length(p_grid));
samples[1:5]
5-element Array{Float64,1}:
 0.715
 0.613
 0.901
 0.49 
 0.692

snippet 3.3

Draw 10000 samples from this posterior distribution

N = 10000
samples = sample(p_grid, Weights(posterior), N);
10000-element Array{Float64,1}:
 0.615
 0.522
 0.872
 0.694
 0.709
 0.851
 0.723
 0.632
 0.84 
 0.778
 ⋮    
 0.361
 0.502
 0.745
 0.456
 0.76 
 0.417
 0.778
 0.349
 0.56 

In StatisticalRethinkingJulia samples will always be stored in an MCMCChains.Chains object.

chn = MCMCChains.Chains(reshape(samples, N, 1, 1), ["toss"]);
Object of type Chains, with data of type 10000×1×1 Array{Float64,3}

Iterations        = 1:10000
Thinning interval = 1
Chains            = 1
Samples per chain = 10000
parameters        = toss

parameters
      Mean    SD   Naive SE  MCSE   ESS 
toss 0.6362 0.1411   0.0014 0.0013 10000

Describe the chain

describe(chn)
Iterations        = 1:10000
Thinning interval = 1
Chains            = 1
Samples per chain = 10000
parameters        = toss

Empirical Posterior Estimates
────────────────────────────────────────
parameters
      Mean    SD   Naive SE  MCSE   ESS
toss 0.6362 0.1411   0.0014 0.0013 10000

Quantiles
────────────────────────────────────────
parameters
      2.5% 25.0%  50.0% 75.0% 97.5%
toss 0.341 0.542 0.6455  0.74 0.882

Plot the chain

plot(chn)
0 2500 5000 7500 10000 0.2 0.4 0.6 0.8 toss Iteration Sample value 0.2 0.4 0.6 0.8 1.0 0.0 0.5 1.0 1.5 2.0 2.5 toss Sample value Density

snippet 3.4

Create a vector to hold the plots so we can later combine them

p = Vector{Plots.Plot{Plots.GRBackend}}(undef, 2)
p[1] = scatter(1:N, samples, markersize = 2, ylim=(0.0, 1.3), lab="Draws")
0 2500 5000 7500 10000 0.0 0.2 0.4 0.6 0.8 1.0 1.2 Draws

snippet 3.5

Analytical calculation

w = 6
n = 9
x = 0:0.01:1
p[2] = density(samples, ylim=(0.0, 5.0), lab="Sample density")
p[2] = plot!( x, pdf.(Beta( w+1 , n-w+1 ) , x ), lab="Conjugate solution")
0.00 0.25 0.50 0.75 1.00 0 1 2 3 4 5 Sample density Conjugate solution

Add quadratic approximation

plot(p..., layout=(1, 2))
0 2500 5000 7500 10000 0.0 0.2 0.4 0.6 0.8 1.0 1.2 Draws 0.00 0.25 0.50 0.75 1.00 0 1 2 3 4 5 Sample density Conjugate solution

End of 03/clip-02-05.jl

This page was generated using Literate.jl.